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Abstract 

A novel representation of functions, called generalized Taylor form, is applied 
to the filtering of white noise processes. It is shown that every Gaussian colored 
noise can be expressed as the output of a set of linear fractional stochastic differ- 
ential equation whose solution is a weighted sum of fractional Brownian motions. 
The exact form of the weighting coefficients is given and it is shown that it is 
related to the fractional moments of the target spectral density of the colored 
noise. 



1 Introduction 

In this paper a novel procedure to represent a stationary Gaussian process by filtering a 
Gaussian white noise process is reported. 

Linear stochastic differential equation excited by a Gaussian white noise process acts as 
a filter, returning an output that is a Gaussian process with Power Spectral Density (PSD), 
or correlation function, related to the equation's parameters. The system linearity, of course, 
implies the normality of the output process. Depending on the structure and on the coefficients 
of the filter equation, processes with different PSD might, at least theoretically, be obtained. 
In practice, such a task is not simple at all, and many papers in literature deal with the 
characterization of filter in order to fit some target PSD. 

Such a spread interest depends on the fact that many real phenomena of engineering and 
physical interest are indeed modeled as stationary Gaussian processes with PSD that is known 
from experimental works. To cite just fews, the most used in the fields of earthquake, wind 
and ocean engineering are the Tajimi-Kanai [20J, the Davenport and the Kaimal [6], [10] and 
the Pierson-Moskowitz |14j spectra, respectively. Representation of such processes as output 
to linear differential equations excited by white noises is a key issue in dynamical analysis 
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of single and multi degree of freedom systems. In earthquake engineering readers can find 
some example in [2], [8], [7] while in the context of offshore design, shaping filters for random 
response analysis are studied in |21j . Auto Regressive {AR), Moving Average {MA) and their 
combination {ARM A) models have been extensively used by many authors to represent colored 
processes. In [TH] the applicability of AR, MA and ARMA algorithms to represent the wave 
motion following a Pierson-Moskowitz spectrum is discussed. Further, in [19] an analog filter 
approximation is presented for the Jonswap spectrum. 

In a spectral fitting problem, a common issue of such spectrum approximations is that the 
determination of filter parameters defining the PSD digital model, relies on some optimization 
criterion. Indeed, once the AR, MA or ARMA model is selected, the coefficients defining the 
pulse transfer function that must fit the target PSD lack of any further interpretation. 

Recently, the first two authors introduced the fractional calculus and the generalized Taylor 
expansion involving Fractional Spectral Moments (FSM) to describe both PSD and correlation 
function in the whole domains ] — oo < oj < oof and ] — oo < r < oo[. That is, the spectral 
fitting problem is easily recast in finding integrals of the target function, which have the 
meaning of complex moments, and represent the coefficients of a series representation. This 
approach based on fractional calculus is of great generality and has also been presented for 
density estimation in probability in [3] and [5]. 

The new issue presented in this paper, extensively based on the latter concepts, is to find 
out the differential equation coupled to the spectral representation based on complex moments. 
In particular, it will be assumed to deal with some spectral data coming from experiments, 
i.e. the target PSD. Then, by the fractional spectral moments two results will be presented. 

The first result is a representation of the stationary Gaussian process. Indeed, it is shown 
that a process with assigned target density can be represented once the FSM are known, by 
means of an expression that involves fractional derivative of a Gaussian noise. Such processes 
have been recently shown to be fractional Brownian motions (fBm) [13]. Then, every Gaus- 
sian colored noise can be thought as superpositions of weighted fBm, and the weights are 
determined by the FSM. 

From this result a further second very remarkable result is achieved. Indeed, by proper 
application of the fractional calculus, a linear fractional stochastic differential equation, whose 
solution is the colored noise process with target density, is found. 



2 A new representation formula for Fourier pair 
functions 



In this section, some preliminary concepts and definitions on fractional operators are summa- 
rized for clarity's sake as well as to introduce appropriate notations. 

Let us consider a Fourier transformable function f{t) and let us denote (f{co) its Fourier 
transform, that is 



Let us recall the definitions of the Riesz fractional integrals and derivatives as follows 



(1) 




and its inverse transform is written as 
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where ^(7) = r(7)cos(77r/2) with T (•) is the Euler gamma function and 7 = p + ii], p > 0, i] G 
R. Their Fourier transforms, in case < p < 1, (|15j. p. 2 17) are 

(4a) ^{(r/)W;w} = M -7 ^{/(*);W 



(4b) T{(V'yf)(t) ] uj} = \uj\^{f(t) ] uj} 

Comparing eqs. (|4ap - (|4b|) it can be simply verified that the relation (D 7 /) (t) = (J -7 '/) (t) 
holds true. Readers should keep in mind that such condition is valid at least for Fourier 
transformable functions here considered. This property allows to calculate the fractional 
integral in eq. (|3ap by simply changing the sign inside the definition of eq.(|3b[) and vice- versa. 
It is to be stressed that this fact is not trivial at all. Indeed, dealing with functions that 
are not Fourier transformable, in general (I> 7 /) (t) 7^ (/~ 7 /) (t) (see [TS], p. 214 and p. 112-3, 
Lemma 5.2). 

For brevity's sake, definitions of the other fractional operators as Riemann-Liouville (RL) 
fractional integral and derivative and of the Marchaud fractional derivative, with their Fourier 
transforms are provided in Appendix Al. Readers can find a complete theory on such operators 
in the excellent monograph of Samko et al. |15j . 

For functions that are Fourier pairs, as those considered in this paper, the eqs.@ are 
very useful to calculate the fractional operators in an easier way with respect to definitions in 
eqs.©. Indeed, by applying inverse Fourier transform to eqs.flD), it leads to 

(5) (Pf) (t) = {V~y)(t) = — J \u>r V M e-^dw 

From the above equation it can be observed that the convolution integrals (|3a|) and (I3bl) 
may be evaluated by making firstly the Fourier transform of f(t), namely ip (uj) and then 
making the Fourier transform of |w| -7 ip (u). This remark is very suitable especially when the 
Fourier transformation of f(t) is known in analytical form. 

Eq.© evaluated in t = assumes the particular meaning of fractional moments of <p(u;) 
in the form 

(6) M (- 7 ) = / 27r(/ 7 /)(0)= f°° | W |">Hd W 

It has been recently shown [3], that the quantities in eq.© are able to represent both 
the function f(t) and (p(co) for symmetric real function. Since the goal of the paper is to 
represent univariate processes as output of filtered white noise and the Gaussian process is 
fully characterized in the probabilistic setting by the (symmetric) PSD, in the following, for 
readability's sake we suppose that f(t) is symmetric. This leads to simplifications, as the 
Riesz integral definition in eq.([3]) can be rewritten in the form 

1 r 00 

(7) (I J f)(0) = —- r 7 ~V(r)dr 

v (7) Jo 

and directly interpreted as Mellin transform (see Appendix A2). In this paper we will deal 
exclusively with f(t) G R and symmetric and consequently (p(co) G R is also symmetric. 
Extension to more general conditions is given in Appendix Al. Then, the corresponding 
representations of the symmetric functions f(t) and (p(uj) are given by 

(8a) f^ = T~ -^V(-7)IC 7 d 7 
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p+ioo 



7-1 



(8b) 



A* (-7) l w 



(I7 



Both the integrals are performed along the imaginary axis with fixed real part p belonging 
to the so-called fundamental strip of the Mellin transform. A deeper insight on Mellin trans- 
form related concepts and more information on the application to fractional derivatives, here 
omitted for readability's sake, are reported in Appendices Al and A2. 

Keeping in mind eq.([6]), eqs.([8]) can be understood as a Taylor integral expansion, because 
by the knowledge of the fractional integrals (or derivatives) in t = 0, it may fully reconstruct 
the function in the whole domain. For this reason, eqs.([SD will be indicated as generalized 
Taylor integral forms. This particular form involving fractional moments of Fourier transform 
of f(t) is entirely new at our best knowledge. It is to be remarked the //(— 7) is able to 
represent both f(t) and its Fourier transform in the whole domains. 

Eq.([S|) assigns further a geometrical and physical meaning to fractional operators in the 
Taylor integral form, relating them with the concept of complex moments of the function 
<p(cu)- Applications of eqs.([8]) for the case of correlation functions and PSD may be found 
in [3], where the fractional moments of the power spectral density assume the meaning of 
fractional spectral moments and were labeled by A(7). 

It is important to observe that in the generalized Taylor expansions ([8j), the integral is 
performed along the imaginary axis and, under the hypothesis that the direct Mellin transform 
in eq.([SD exists, such integrals do not diverge in virtue of the Mellin inverse theorem. Moreover 
it has to be stressed that the integration does not depend on the particular choice of p except 
for the limitation, as previously outlined, that p shall belong to the fundamental strip of the 
Mellin transform. This may be explained by the fact that the integrand in the fundamental 
strip is holomorphic. More information on this topic may be found in [3]. 

Hereinafter, taking full advantage of the results presented in this section, the fractional 
filter representing a given target PSD will be found. 



3 Transfer function representation by H-Fractional 
Spectral Moments 



Objective of this paper is to represent a normal stationary process with assigned power spectral 
density as the output of a fractional differential equation. To this aim, let us consider a 
Gaussian white noise process, W(t) with zero mean and correlation function E[W(t)W(t + 
t)] = q 8(t), and power spectral density Sw = 9/(27r), where q is the intensity parameter. Let 
us indicate the ideal linear system as 



where C (•) is a linear differential operator applied to the response Y(t). The solution to 
eq.Q can be characterized by the impulse response function h{t) and its Fourier transform 
H(cj) namely the transfer function, and can be expressed by Duhamel integral 



Indicating by 5y(w) and Sw(w) the power spectral density of the output and the input, 
respectively, from the linearity of the system it follows 



where Sy(oj) is the target spectrum. Eq.Q may be considered as a filter of the white noise 
process. Let us now suppose that the differential operator C (•) is such that Sy(w) overlaps 



(9) 



C(Y(t)) = W(t) 



(10) 




(11) 



S Y (w) = \H M| 2 S w (w) = ^-\h H| 2 
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an assigned PSD arising from a physical phenomenon. In order to find the unknown operator 
C (•) we assume Arg {H (w)) = 0, that is, from eq. (jll[) the transfer function can be expressed 
as 

(12) H{ U ) = \H{u)\ = J—Sy(u) 

Correspondingly, as the transfer function is assumed to be real, the impulse response 
function of the linear system C(Y(t)) remains symmetric. That is, by enforcing the condition 
Im {H (uj)) = we get a non causal differential equation. In spite of the causality condition is 
violated, the output of eq.([U]) remains a strictly stationary Gaussian process. In the following, 
we will firstly introduce the fractional moments of the function H(u); then we will represent 
H(uj) as sum of fractional moments and finally the expression of the process Y(t) with target 
PSD Sy(w) will be given in terms of the fractional moments of H{oo). As first step, in order 
to represent the transfer function H(oj), let us define the fractional moment of H{oj) labeled 
in the following as Uh ( — 7) £ C, that in virtue of eq.([6]) is written as 

(13) U H (-7) d = H (w) doj, Re 7 > 

that will be called H-Fractional Spectral Moment (H-FSM) function. By means of eq.®, the 
H-FSM function allows us to fully reconstruct both the function H(uS) and h(t) in the whole 
domains. The H-FSM are fractional integrals and derivatives of h(t) evaluated in zero (as 
given in eq.dSP), that is 

(14) 2tt (.PTi) (0) = 27r(P~ 7 /i)(0) = n H (-7) , Re 7 > 

Specifying the eq.([5]). previously written for two general Fourier pair, we obtain the rep- 
resentation of the impulse response in the time domain as 

^ r-p+ico 

(15) h(t) = —^ j/(7)lM-7)i~ 7 d 7 , t>0, < Re 7 < 1 

(2tt) l Jp~ioo 

and the transfer function in the form 

^ rp+ico 

(16) H(oj) = — IL H (-7)|w| 7_1 d7, 0<Re 7 <l 

J p—ioo 

It might be useful to consider the H-FSM as a third representation in the 7 variable of the 
response of the dynamical system under study. 

The integrals involved in (|15p and (|16p can be approximated in discrete form, operating 
a truncation in the r\ axis. Indeed, calculating the integral up to a certain value, called fj, 
dividing [—fj,fj] in 2m intervals of amplitude A?7 = fj/m, with m G N and evaluating the 
integrals at the values Jk = P + ikArj, eqs. ()15p and ()16p can be approximated in the form 

A m 

(17) h (t) - — ^ Yl v (Tfc) n ^ ("7fc) t~ Jk 

(2tt) , m 



(18) H (p)* Y, Uh M 7 ' 5 " 1 

k=—m 

A wider discussion on the truncation of the integrals performed along the imaginary axis 
may be found in [3]. 

Now, having represented the transfer function both in exact and approximated form in 
terms of H-FSM we are ready to infer an analytic expression for the process Y(t) with target 
PSD. 
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The input-output relation for linear system in Fourier domain is written as Y (ui, T) = 
H (lu)W (lu,T), where T > is a truncation bound. Of course this relies on the fact that 
W(t) is a stationary noise. Thus, bearing in mind eq. (|16j) one obtains 

i r-p+ioo 

(19) Y(oj,T) = — U H (-7) M 7 " 1 W (u/, T) d 7 

J p—ioo 

due to the linearity of the operator C (•). Applying an inverse Fourier transform, the response 
to the linear system assumes the form 

i rp+ioo . . 

(20) Y (t) = — U H (- 7 )jF- 1 \\u>r 1 W(cj,T);t\d 7 

47T2 J p-ioo < > 

By recalling eq.gaD, it follows that lim J 7 " 1 (|a;| 7_1 TV (u>, T) ; t\ = (i^W) (t), that 
introduced in the latter, gives 

i rp+ioo 

(21) Y {t) = — U H (-7) (P-IW) (t) d 7 

J p—ioo 

with p > as already found by means of the impulse response function, that is the new exact 
representation of the stationary process with assigned PSD. 

Some comments are necessary to highlight the peculiar aspects of this new representation 
of a stationary process, i) The resulting process Y(t) is reconstructed by the knowledge of the 
H-FSM previously defined in eqs.(|13p. calculated on the transfer function given by eq. (|12p . 
through the target power spectral density Sy(uj). There is no need of using optimization 
criteria like in ARM A, or similar models, because the coefficients figuring in the representation 
have the meaning of being complex moments of the H(u). ii) Eq. (|21|) is an integral along 
the imaginary axis with fixed real part p, belonging to the interval [0,1]. The choice of p 
inside this interval does not influence the integral because the integrand is holomorphic inside 
such a interval. Hi) Fractional integrals of white noise processes have recently attracted many 
authors. In fact, such operators are connected with the so called Fractional Brownian noises 
and motions as reported in [1] , [9] . Although white noise processes have nowhere differentiable 
path, the operation of fractional integration and derivation is indeed meaningful. To get a clue 
on this, it suffices to recall the definitions in eqs. ([3a]) - (|3bp which stress the convolution nature 
of the fractional operators. The kernel of the integral smooths the singularity of the process 
path in a such a way that it is therefore well-defined. In ([22], pp. 65- 70) very interesting 
results on a class of nowhere differential function (the Weierstrass function) is also presented. 
iv) The process Y(t) as expressed in eq. ([2"Tj) is suited to be computed by proper discretization 
of the integral involved. Indeed, truncating the integral up to a certain value, called fj, and 
dividing [—77,77] in 2m intervals of amplitude Ar/ = fj/m, the integral can be approximated by 
the sum 

A m 

(22) Y (t) = - Jf £ U H (- 7fc ) (I'-^W) (t) 

k=—m 

with 7fc = p+ikArj (0 < p < 1). This approximation carries out a truncation and discretization 
error, that can be made arbitrarily small. Eq. (|2ip . or its discretized counterpart given in 
eq. (|22p . are the new representation of the process Y(t) whose PSD matches the target one. In 
particular eq. ([22"j) shows that the Y(t) may be obtained as the superposition of the fractional 
integrals of the white noise process. 

In order to validate eq.(|22p. we will show that with some simple steps it coincides with 
the well-known Shinozuka's representation [16]. Let us consider a band limited white noise 
process with one-side PSD 



(23) G w H 



q/ir < uj < oj 
otherwise 
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where oj is some cut-off frequency and let be Aoj a discretization of the oj axis such that 
Aw = oj/n, then the spectral representation theorem [16] 



(24) W (t) = Jim V^oN sin fat + fa ) 

holds true, where fa are realizations of a random variable uniformly distributed in [0, 2ir]. The 
fractional integral ( y I 1 ~ lk W) (t) can be directly computed by Mathematica and is 

n 

(25) (I 1 " 7 * W) (t) = lim^ V 2A ^ 9 A" N 7fc_1 sin fat + fa) 

Aw->0 j=l 

Introduction of the latter in eq. ([2Tj) leads to 

(26) Y (t) = lim £ ^A^/vrsin fat + fa) — / n H (-7) fa\~< " x d 7 

that, bearing in mind eqs. (|12p and f)16[) . is simplified in the form 



(27) F (t) = Jim Y, V 2G V Aw sin M + 0i 



.. ->oo _ 
Ao;-»0 j=l 



where Gyfa = 2SyfaUfa is the target unilateral power spectral density. Eq. ([27|) is the 
well known Shinozuka's representation |16| for a process with assigned spectrum, and then 
the process reconstruction by eq. ([22|) is proved. 

3.1 Two Relevant Examples 
3.1.1 Pierson-Moskowitz spectrum 

Consider the unilateral Pierson-Moskowitz (PM) spectrum [T5] Gpm given by 



(28) G PM fa 



^exp(-^); oj>0 
0; oj < 



with assigned parameter c. Let us study the linear system excited by a Gaussian white noise 
process as indicated in eq.© such that Y(t) has a power spectral density with the assigned 
form Sy(±o;) = ^Gpm(w), with oj > 0. It follows from the definition in eq. (|13|) that the 
H-FSM can be easily evaluated by Mathematica and assume the form 



(29) n H (-7) = 2t + ^5-§-ic 5 / 8 v / ^/-r 



3 7 
8 4 



having introduced HpMfa = a/vt \Gpm fa)\ /<Z, as given by eq. (fTT|) . 

The approximated impulse response function hpM(t) an d the transfer function H(oj) of 
the system under exam have been found by applying eqs. (fTT|) and (fT8|) and plotted in Fig.© 
and Fig. ([2]), respectively. The parameters used are: c = 2.72, p = 0.5, m = 25, Ar] = 0.6 
(fj = 15). It can be noted that the comparison between the exact functions (continuous line) 
and the approximated one (dotted line) is very good. As confirmed by the log- log plot in 
Fig-©i the two functions coincide in a very wide interval. Although in Fig.© the impulse 
response function is plotted only for t > 0, it is a symmetric function because the differential 
equation associated to eq. ([T7|) is non causal. 
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Figure 2: Approximate Hpm(oj) (dotted) contrasted with the exact one 



3.1.2 Davenport's spectrum 

In wind engineering, the Davenport's spectrum 

A7rk V r 2 q D (oj) 2 



(30) S D (u) 



0J\ 



4/3 



is widely used to represent wind fluctuation, where V r is the mean wind speed at the reference 
level, ko is a roughness characteristic of the analyzed site and <7_d(^0 = 1200aj/(27rT4.). For this 
application we selected ko = 0.01, Vr = 15m/s, q = 1 as spectrum parameters and p = 0.7, 
m = 30, Ar] = 0.5 {f) = 15). In fig.(jl|) the approximated transfer function Hr>(uj) is contrasted 
with the exact one defined from eq. (|30p . With this further example we want to stress that the 
spectrum behavior in the neighborhood of the zero has no influence on the applicability and 
the efficiency of the method. Indeed both the flat and the steep behavior of the PM and the 
Davenport spectra, are very good approximated by H-FSM. 



4 Fractional differential equations of the linear filter 

The method presented in the previous section characterizes the process Y(t) with assigned 
spectrum, by the exact expression in eq. ([2"T]) or by the approximated one in eq. ([2"2"jh In this 
section, the differential equation whose response is Y(t) will be found in approximated form, 
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Figure 3: Log-Log plot of the approximate form of /ipif(t) (dotted) contrasted with 
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Figure 4: Approximate H D (u) (dotted) contrasted with the exact one 



making extensive use of composition properties of fractional operators reported in Appendix 

m 

A3. Firstly, let us rewrite eq. ([22|) in the form Y (t) = Yl ^fc (*) where 

k=—m 

(31) Y k (t) = a k (l 1 -^W)(t) 

with Ofc = ArjH}{ (— 7&) / (47r). Applying to both sides of the latter equation the operator 
D 1 ~ Jk and exploiting the composition rule reported in Appendix A3 leading to eq. fTTj) . we 
obtain 

(32) (D^Yk) (t) = a k W (t) 

Each component of Y(t) is therefore the solution of a linear fractional stochastic differential 
equation excited by a Gaussian white noise. As the noise is the same for each Y k (t), the k 
components are dependent each other. For a complete description of the whole process Y(t), it 
suffices to perform a time derivative of the first order to eq. (|31|) and summing for k = —m, m 
obtaining 

c\Y m 

(33) — = £ a k (D^W)(t) 

k=—m 

that is the linear fractional stochastic differential equation searched. 
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In order to provide the probabilistic characterization of the response process Y(t) let us 
consider the truncated Fourier transform of eq,(|32p yielding 



1-7* 



(34) 

or in terms of PSD 

(35) S Yk (u;) = \a k \ 2 \lo\ 2 ^^ S' 
while the cross spectral density takes the form 

(36) S Yh Yj (w 



Y k (u>,T) = a k W(u,T) 



I 2 I i2(7 fe -l) 
\ a k\ M 



(q/2n) 



* i i yic — l m i 7t — l 
g^g^- |cj| /ft ( I a; i 



(?/2tt) 



where the symbol * means complex conjugate. 

It follows from eq. ()36p that the PSD of the target process Y{t) is written as 



(37) 



Sy (u)= Yl SY k Y i H 



-m j=—m 



Therefore, the cross-spectral density between the dependent components Y k characterizes 
the target power spectral density Sy(u;). We want to stress that each component Sy k gives 
a contribute in defining the total power spectral density Sy(u;) and that only the whole 
ensemble returns a very good approximation. This fact is highlighted in Fig. ([5]) where the 
approximation of a Pierson-Moskowitz spectral density, Sy(uj), by components' cross-spectral 

m r 

densities Sy (w) = Yl Yl ^YkYj varying r is plotted. This figure shows in particular 

k=—m j=—m 

that every single component Y k (t) in eq. (|3ip gives a contribution, because the convergence is 
attained only when r = m. Summing up, every component Y k weights in the reconstruction 
of the target PSD. 
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Figure 5: Approximation of the Pierson-Moskowitz spectrum Sy(u), by components' 
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5 Comparison with time series models 



The processes Y(t) and W(t) up to now considered are continuous in time. In this section we 
show how to express the process Y(t) with target spectral density in a sampled time space. 
To this aim, let us discretize the time axis in equally spaced intervals of amplitude At > 0. 
With the notation Y ktt _j, j G N we denote Yf. (t — jAt) and consequently Y k; t is the process at 
the instant of time t. The same notation applies to the noise such that Wt-j = W (t — jAt). 

The discretization of the Riesz operator can be tackled by the Griinwald-Letnikov approach 
as reported in [15] that reads 



(38) Df>Y k)t = hm [ Y *j (P) Y Kt-j+ Y X J (Z 3 ) Y W 

=o j=0 



with 

(-iy 

(39) Xj (/?) 



2AtP cos (/Stt/2) 

Eq. (|38|) is very useful in numerical applications because for a small value of At > gives 
an approximation formula of the first order. For our scope eq. (|38p is extremely useful because 
gives the possibility to interpret the fractional filter equation in eq. (|32p also in the discrete 
time domain in order to make a fruitful comparison with other well-known time series model. 

Then, by means of the discretization in step of small At > and substituting the approx- 
imated Riesz operator of eq. ()38p in eq. (|32p . the following time series is obtained 

oo 

(40) Y A i ( X " T*) { y M-i + Y k,t+j} = a k W t 

3=0 

This representation highlights that this fractional time series model is determined by the 
knowledge of the past of the process in the first term and the future of the process in the 
second term. This happens because the filter is non-causal. From eq. (|40p it is possible to 
find the transfer function by applying the z-transform. Remanding to the textbook of [12] on 
details concerning the z-transform, it can be proved that the pulse transfer function for the 
single component Yfc is 

(41) H k (z) 



j=0 



The fractional filter characterized by eq. (|40p and (j41|) is suitable to be compared with the 
classical Auto Regressive model of order p £ N, denoted as AR(p), whose time series has the 
form 

p 

(42) y t + Y a jyt-j = GW t 

i=i 

and pulse transfer function given by the equation 

(43) H AR{p) (z) = ^ 

1 + £ a 3 z~i 
j'=i 

Comparing the pulse transfer function in eq. (|4ip and in eq.(|43p it is worth to note that 
both expressions present the z variable only in the denominator, sharing the same algebraic 
structure. Remarkable difference is that the coefficients in the fractional filter are evaluated 
by the H-FSM in exact form while in the AR model a Jule- Walker scheme must be adapted 
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for finding the coefficients a^; by contrast once the H-FSM are evaluated, all the coefficients 
in the fractional differential equation are readily found. 

Same kind of analogy can be pursued starting from eq,(|31|) that assumes the form 

oo 

(44) Y k)t = a k J2 (-0) {W t -j + W t+j } 

3=0 

in the discrete time model and with transfer function 

(45) H k (z) = a k Y, Aj (-0) { z~ j + z J } 

3=0 

Indeed it is easy to recognize that eq. (|45|) is conceptually similar to a Moving Average model 
of order q € N indicated by MA(q), characterized by the temporal series 

g 

(46) y t = J2 b jW t - j 

3=0 

and transfer function 

g 

(47) H MA(q) (z)=^2b jZ ^ 

3=0 

Comparison between MA model and digital filter obtained by the proposed procedure also 
in this case reveals the same structure. It is to be stressed that for MA model, the coefficients 
bj are evaluated by some optimization criteria while in eq. ()45|) the coefficients a k are evaluated 
by H-FSM. 

The last step we develop regards the ARM A model and follows plainly. First of all, let 
us exploit the composition rule given in Appendix in eq. (f54bl) into eq.([32]), applying D"' 3 " 1 to 
both sides of eq. (|32p . In this way the filter equation reads 

(48) (£>T'-T*y) (t) = ak {pf'^W) (t) 
that in the discrete time form is 

oo oo 

(49) X J (7- " 7*) {Yk,t-j + Y kjt+j } = a k Y, Aj (1 " 7.) {Wt-j + W t+j } 

3=0 j=0 

and whose transfer function is 

oo 

(50) H k (z) = 

E Aj (7, - 7fc) {z-i + zi} 

3=0 

The latter should be compared with the ARM A model, that is characterized by the time 
series model 

(51) X>Y^=E C ^-; 

3=0 j=0 

and transfer function 

g I v 

(52) Harm A = z ^ / d i z ~' 

3=0 I 3=0 

Also in this case direct comparison between eq. ([50|) and eq. ([52|) shows the correspondence 
between the two representations. 
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6 Conclusions 



This paper sheds a new light in the representation of stationary Gaussian colored noises. It 
is shown that every stationary Gaussian process with assigned power spectral density is equal 
to an integral involving the H-Fractional Spectral Moments (H-FSM) and Riesz Fractional 
integrals of the Gaussian white noise process. In this sense, the first result of this paper is a 
new exact representation for stationary processes. 

It has been further shown that the H-FSM can be easily evaluated, and two pathological 
applications like the Pierson-Moskowitz and the Davenport density spectra have been reported. 

Moreover, properly approximating the exact integral representation in to a series form, a 
set of fractional stochastic differential equations excited by an external white noise process 
has been provided. Its solution has the desired target power spectral density. It has also 
been shown that every stationary Gaussian process is approximated by a fractional stochastic 
differential equations, that is a stochastic differential equation excited by fractional derivatives 
of a Gaussian white noise. 

The representation proposed has been rewritten in time series form and compared with 
AR, MA and ARM A time series. 

A Appendix 

The generalized Taylor form applied in the paper to symmetric functions expressed in terms of 
Riesz operators is valid also once this condition is dropped, as reported in [3]. For readability's 
sake we report the main results, introducing the proper operator of Riemann-Liouville and 
Marchaud. 

A.l Riemann-Liouville (RL) and Marchaud fractional opera- 
tors 

We recall the definitions of the Riemann-Liouville (RL) fractional integral and derivative, 
(l±f) (t) and CD±f) (t), respectively, and of the Marchaud fractional derivative, (D±f) (*)> 
given by 

(53a) (4/) (t) = f J-j J™ C'-'f (t T d£ 

(53b) M/)W ^_^L_d^^ /(tT0 ^ 

(53c) (did (t) = f f i-j J™ r 7 " 1 (/ (* t o - / it)) de 

where T (•) is the Euler gamma function and 7 = p+ii], p > 0, r] £ R. The Marchaud definition 
of the fractional derivatives is more convenient with respect to the Riemann-Liouville in the 
sense that they exist also for function that do not vanish at infinity, and for this reason it will 
be preferred. 

Fourier transforms of the fractional derivatives and integrals above defined are 
(54a) T { (4/) (t) ; u} = (t^H ^ {/(<); "} 

(54b) T { (Dlf) (t);oj} = ( T iu^ T {/ (t) ; u] 
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Applying same reasoning reported in section 1 for the Riesz operators, comparing eqs, (|54ap - 
(|54bp it can be simply verified that for functions that are Fourier transformable (D±f) (*) = 

fj (t), although not valid in general. Then, by making an inverse Fourier transform of 

the eq. (|54p the following representations 

(55) (4/) (t) = (D^f)(t) = — J ( Tiw )-T ,p ( W ) e-^doj 

holds true. Moreover, eq . (|55|) evaluated in t = assumes the particular meaning of fractional 
moments of <p{ui) in the form 

(56) p T (-7) = 2n (4/) (0) = f°° (=Rw)~ 7 <f («,) du 



Once the fractional moments have been introduced, let us interpret the fractional operators 
in eq. ()56p as Mellin transform of the function f(t) 

/•oo 

(57) (4/)(o)r( 7 )= / r-VCTOtf 

J 

as suggested in [15]. Then, making an inverse Mellin transform, the following representation 
of the function f(t) 

i i-p+ico 

(58) / (Tt) = — r ( 7 ) (4/) (0) *" 7 d7 

holds true, where the integral is performed along the imaginary axis with fixed real part p 
where p belongs to the fundamental strip of the Mellin transform. 

By inserting eq. (|56p in eq. (|58|) . the integral representation for both f(t) and its Fourier 
transform may be written as 



i fp+ ioc r fa) 

(59a) / ( T t) = — / -f>-p T (- 7 ) i-^7 

i/ p — 200 



p+ioo 

d7 



(5%v M = 7TTTT / r ( 7 ) r (7 - 1) M _ (-7) M 7 " 1 +^+ (-7) (-H 7 " 1 



(2vr) z i 

A. 2 Some concepts on Mellin transform 

The Mellin transform ([T3], p.25; |17j . p. 41) of a function / (x), x > 0, s £ C is defined as 

/>oo 

(60) ^ (5) = M {/ (x);s}= f (0 f^dt 

Jo 

along the inverse operator 

PC+IOD 

(61) f{x) = M~ 1 {ip{s);x} = {2wi)' X ip(s) X - s ds 



with c = Re (s) 

From the definition, it follows that the convergence of the Mellin transform depends on 
the behavior of the function f(x) at zero and infinity, that is, assuming 



(62) f{x) 



O (a?P) , x -S- 
O (x q ) , x — > 00 



the complex function (/? (s) is analytic inside the so called fundamental strip —p < Res < —q. 
The inverse Mellin transform must be calculated selecting c = Res inside the fundamental 
strip. Tables of Mellin transforms of commonly used functions are given in ([IT], p. 527), while 
for a complete theory readers are referred to [TT]. 
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A. 3 Compositions rules for Riesz fractional derivatives 

By manipulating the definition of the Riesz fractional operators, it can be proved that the 
following identities 

(6to) (^)W = 5=^{(^)W + (^)W} 

(<Bb) (° 8/ ) <«> = Y^wm i K') w + ( D - { ) w} 

in terms of RL fractional integrals and Marchaud derivatives hold true. 

The left handed \I+f) (t) and the right handed [I-fj (t) operators, are connected by 



the equations 

(64a) (itf) (t) = cos OStt) (/£/) (t) + sin (/3tt) ft { (/J/) (s) ; t} 

(64b) (ij/) (t) = cos GStt) (t) - sin (/3vr) ft { (/!/) ( S ) ; t} 

where ft is the Hilbert transform defined as 

(65) n{f(s);t} = - [°° f^ds 

The eqs. (|64p are useful to find the so called composition rules of fractional operators, 
which allow to simplify forms of the type I5"j£/. In particular, applying and D^_, to 
eq. (|64ap and eq . (|64b|> . respectively, and taking into account 

(66) (Di4f) (<) = / (*) (Dtlif) (t) = f (t) 



it follows 



(67a) W (t) = cos (/3vr) / (t) + sin (/3vr) ft {/ (s) ; *} 



(67b) (D^_4f) 
having used the property 

(68) H(I p ± f)(t) = [l p ± Hf)(t) 



Composition rules involving the Riesz operators can now be worked out based on the 
previous properties, showing that D"l^ f = f. Indeed, from the definitions in eq. (|63D 

(69) D^/ = ^ 1 (D^lif + Btltf + D%ff + Dtlif 

(2 cos (p7r/2)J v 



that using eq. (|66|) , (jSTj) and 

(70) 2 + 2cos W = 1 

(2cos(/3vr/2)) 2 

gives the result searched 

(71) D^Pf = f 

Other compositions rules and applicability criteria of the formulas in this appendix are 
adapted from ([E], Chapter 3). 
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